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Abstract 

The aim of this tutorial survey is to revisit the basic theory of relax- 
ation processes governed by linear differential equations of fractional 
order. The fractional derivatives are intended both in the Rieamann- 
Liouville sense and in the Caputo sense. After giving a necessary out- 
line of the classical theory of linear viscoelasticity, we contrast these 
two types of fractional derivatives in their ability to take into account 
initial conditions in the constitutive equations of fractional order. We 
also provide historical notes on the origins of the Caputo derivative 
and on the use of fractional calculus in viscoelasticity. 
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Introduction 



In recent decades the field of fractional calculus has attracted interest of 
researchers in several areas including mathematics, physics, chemistry, engi- 
neering and even finance and social sciences. In this survey we revisit the 
fundamentals of fractional calculus in the framework of the most simple timc- 
depcndcnt processes like those concerning relaxation phenomena. Wc devote 
particular attention to the technique of Laplace transforms for treating the 
operators of differentiation of non integer order (the term "fractional" is 
kept only for historical reasons), nowadays known as Riemann-Liouville and 
Caputo derivatives. We shall point out the fundamental role of the Mittag- 
Lcfflcr function (the Queen function of the fractional calculus), whose main 
properties are reported in an ad hoc Appendix. The topics discussed here 
will be: (i) essentials of fractional calculus with basic formulas for Laplace 
transforms (Section 1); (ii) relaxation type differential equations of fractional 
order (Section 2); (iii) constitutive equations of fractional order in viscoelas- 
ticity (Section 3). The last topic is treated in detail since, as a matter of 
fact, the linear theory of viscoelasticity is the field where we find the most 
extensive applications of fractional calculus already since a long time, even 
if often only in an implicit way. Finally, we devote Section 4 to historical 
notes concerning the origins of the CapTito derivative and the use of fractional 
calculus in viscoelasticity in the past century. 



1 Definitions and properties 

For a sufficiently well-behaved function f{t) (with t G R^) we may define the 
derivative of a positive non-integer order in two different senses, that we refer 
here as to Riemann-Liouville (R-L) derivative and Caputo (C) derivative, 
respectively. 

Both derivatives are related to the so-called Riemann-Liouville fractional 
integral. For any a > this fractional integral is defined as 

J? m := fit - tT-' fir) dr , (1.1) 

1 [a) JO 

roo 

where r(a) :— / e~" u°'~^ du denotes the Gamma function. For existence 
Jo 

of the integral (1) it is sufficient that the function fit) is locally integrable 
in R+ and for t ^ behaves like Oit~^) with a number u < a. 
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For completion we define J^ — I (Identity operator) . 
We recall the semigroup property 

Furthermore we note that for a > 

■^"^"^ rr^^Il!^ 7>-l- 
r(7 + 1 + a) 



[1.2) 



(1.3) 



The fractional derivative of order > in the Riemann-Liouville sense is 
defined as the operator which is the left inverse of the Riemann-Liouville 
integral of order // (in analogy with the ordinary derivative), that is 



[lA) 



If m denotes the positive integer such that m — 1 < /j, < m , we recognize 
from Eqs. (1.2) and (1.4): 



DU{t):^D^jr'f{t), m-l<ii<m. 
In fact, using the semigroup property (1.2), we have 



(1.5) 



Jit = A"" JT~^ Jt = DT jr = i. 



't 



Thus (1.5) imphes 



' f{r)dT 
V{m - Jo {t - t)/^+i- 

fit), 



, m — 1 < iJ, < m; 
II — m. 



;i.5') 



For completion we define = I . 

On the other hand, the fractional derivative of order /j, in the Caputo sense 
is defined as the operator such that 



.D^ fit) fit) , m - 1< < 



m . 



(1.6) 



This implies 
*Dmt) = { 



1 



r(m - /i) Jo {t - r)^+i 
d'^ 



— , m — 1 < a < m; 



:i.6') 



jj, = m . 



Thus, when the order is not integer the two fractional derivatives differ in 
that the standard derivative of order m does not generally commute with 
the fractional integral. Of course the Caputo derivative (1.6') needs higher 
regularity conditions of f{t) than the Riemann-Liouville derivative (1.5'). 

We point out that the Caputo fractional derivative satisfies the relevant prop- 
erty of being zero when applied to a constant, and, in general, to any power 
function of non-negative integer degree less than m , if its order /i is such 
that m — 1 < < m . Furthermore we note for /i > 0: 



r(7 + i) 
r(7 + 1 - fi) 



7 > -1 



:i.7) 



It is instructive to compare Eqs. (1.3), (1.7). 

In [27] we have shown the essential relationships between the two fractional 
derivatives for the same non-integer order 



m—l 



fit) - E /^'Ho^ 



fc=0 



1 /W(o+)t'=-A' 



m — 1 < fi < m . 



:i.8) 



In particular we have from (1.6') and (1.8) 



1 



r(l - /i) Jo {t - r)^ 



dr 



fit) - /(o- 



A^7(t)-/(o- 



0</i<l. (1.9) 



r(i-/i) 



The Caputo fractional derivative represents a sort of regularization in the 
time origin for the Riemann-Liouville fractional derivative. We note that 
for its existence all the limiting values f^^\0^) := lira /(t) are required 

to be finite for A; = 0, 1, 2, . . . , m — 1. In the special case /'■'^^(O^) = for 
/c = 0, 1, m — 1, the two fractional derivatives coincide. 

We observe the different behaviour of the two fractional derivatives at the end 
points of the interval (m—l, m) namely when the order is any positive integer, 
as it can be noted from their definitions (1.5), (1.6). In fact, whereas for 



m 



both derivatives reduce to D™, as stated in Eqs. (1.5'), (1.6'), due 
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to the fact that the operator = I commutes with D]^, for fi ^ [m — l)"*" 
we have 

r D^fit) ^ nr J] fit) = D^-'^ fit) = f^-\t), 

/X (m - 1)+ : <^ (1.10) 

[ Mfit) ^ JlD^fit) = -/(™-i)(0+). 

As a consequence, roughly speaking, we can say that is, with respect to 
its order /i , an operator continuous at any positive integer, whereas ^D^ is 
an operator only left-continuous. 

The above behaviours have induced us to keep for the Riemann-Liouville 
derivative the same symbolic notation as for the standard derivative of integer 
order, while for the Caputo derivative to decorate the corresponding symbol 
with subscript *. 

We also note, with m — 1 < < m , and Cj arbitrary constants, 

m 

Df/W = Digit) ^ fit)=git) + Y.c,t^-\ (1.11) 

i=i 

m 

Mfit) = .Digit) ^ fit) = git) + Y.c,t^-^ . (1.12) 

i=i 

Furthermore, we observe that in case of a non-integer order for both fractional 
derivatives the semigroup property (of the standard derivative for integer 
order) does not hold for both fractional derivatives when the order is not 
integer. 

We point out the major utility of the Caputo fractional derivative in treating 
initial- value problems for physical and engineering applications where initial 
conditions are usually expressed in terms of integer-order derivatives. This 
can be easily seen using the Laplace transformatioE0 



^The Laplace transform of a well-behaved function f{t) is defined as 

f{s)=C{f{t);s}:= e-'^fit)dt, 3? (s) > a/ . 
Jo 

We recall that under suitable conditions the Laplace transform of the m-derivative of f{t) 
is given by 



m— 1 

£ {or f{t);s} = lis) - ^ s—i-^- /('^)(0+) , /(^■)(0+) lim D^f{t) . 



k=0 
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For the Caputo derivative of order /i with m — 1 < jj, < m we have 

m—1 

C {M f{ty,s} = f{s) - Y: s^-'-' f^'\0^) , 

(1.13) 

/W(0+) := limD^fit). 
The corresponding rule for the Riemann-Liouville derivative of order fi is 

m—1 

C {A" f{ty,s} = m - Y: s--'-' g^'\0^) , 

'^=° (1.14) 

g(>'\0+) := \im D^g{t) , where g{t) := J^"^^ f{t) . 

Thus it is more cumbersome to use the rule (1.14) than (1.13). The rule 
(1.14) requires initial values concerning an extra function g{t) related to the 
given f{t) through a fractional integral. However, when all the limiting values 
y(fc)(^0+) for = 0, 1, 2, . . . are finite and the order is not integer, we can prove 
that the corresponding g^''\0~^) vanish so that formula (1.14) simplifies into 

£{Dr/(t);s} = s^/(s), m-l<fi<m. (1.15) 

For this proof it is sufficient to apply the Laplace transform to the second 
equation (8), by recalling that £{)f:°';s} = T{a + l)/s"+^ for a > —1, and 
then to compare (1.13) with (1.14). 

For fractional differentiation on the positive semi-axis we recall another def- 
inition for the fractional derivative recently introduced by Hilfer, see [67] 
and [115] , which interpolates the previous definitions (1.5) and (1.6). Like 
the two derivatives previously discussed, it is related to a Riemann-Liouville 
integral. In our notation it reads 

^M,- _ ^.(i-M) ^1 ^ 0</i<l, 0<z/<l. (1.16) 

We can refer it to as the Hilfer (H) fractional derivative of order fi and type 
u. The Riemann-Liouville derivative corresponds to the type i/ = whereas 
that Caputo derivative to the type u = 1. 

We have here not discussed the Beyer-Kempfie approach investigated and 
used in several papers by Beyer and Kempfie et al.: this approach is appro- 
priate for causal processes not starting at a finite instant of time, see e.g. 
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[HI EH] . They define the time- fractional derivative on the whole real line as a 
pseudo-differential operator via its Fourier symbol. The interested reader is 
referred to the above mentioned papers and references therein. 

For further reading on the theory and applications of fractional integrals and 
derivatives (more generally of fractional calculus) we may recommend e.g. 
our CISM Lecture Notes [551 EH EO] , the review papers [851 ESI 1116] , and the 
books [HHl [701 [661 EH EZl m (M EM EM with references therein. 



2 Relaxation equations of fractional order 

The different roles played by the R-L and C derivatives and by the inter- 
mediate H derivative are clear when one wants to consider the correspond- 
ing fractional generalization of the first-order differential equation governing 
the phenomenon of (exponential) relaxation. Recalling (in non-dimensional 
units) the initial value problem 

— = -u{t), t>0, with m(0+) = 1 (2.1) 

(Jjv 

whose solution is 

u{t) = exp(-t) , (2.2) 

the following three alternatives with respect to the R-L and C fractional 
derivatives with /i G (0, 1) are offered in the literature: 

^D^ u{t) = -u{t) , t>0, with u(0+) = l. (2.3a) 

u{t) = -u{t) , t > , with lim J^'^ u{t) = 1 , (2.36) 

^ = -DI'^ u{t) , t>0, with u(0+) = 1 , (2.3c) 



In analogy to the standard problem (2.1) we solve these three problems with 
the Laplace transform technique, using respectively the rules (1.13), (1.14) 
and (1.15). The problems (a) and (c) are equivalent since the Laplace trans- 
form of the solution in both cases comes out as 



+ 1 
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whereas in the case (b) we get 



1 s'^^^ 
u{s) = r = l-s -. (2.5) 

The Laplace transforms in (2.4)-(2.5) can be expressed in terms of functions 
of Mittag-LefHer type, for which we provide essential information in Ap- 
pendix. In fact, in virtue of the Laplace transform pairs, that here we report 
from Eqs. (A.4) and (A.8), 

C{E,{-Xn,s} = ^^, (2.6) 

C{t''''E,A~^ns} = ^^, (2.7) 
with /i, G and A G R, we have: in the cases (a) and (c), 

M(t) = ^(t) := E^(-t'^) , t>0, 0</i<l, (2.8) 
and in the case (b), using the identity (A. 10), 

uit) = ^t) := t-^'-^^E,,, i-tn = -j^E, i-t^ , t > 0, < /i < 1. (2.9) 

It is evident that for ^ ^ 1~ the solutions of the three initial value problems 
(2.3) reduce to the standard exponential function (2.8). 

We note that the case (b) is of little interest from a physical view point since 
the corresponding solution (2.9) is infinite in the time-origin. 

We recall that former plots of the Mittag-Leffler function "^{t) are found 
(presumably for the first time in the literature) in the 1971 papers by Ca- 
puto and Mainardi [2 80 and [21] in the framework of fractional relaxation 
for viscoelastic media, in times when such function was almost unknown. 
Recent numerical treatments of the Mittag-LefHer functions have been pro- 
vided by Gorenflo, Loutschko and Luchko [56j with MATHEMATICA, and 
by Podlubny [97] with MATLAB. 



•^Ed. Note: This paper has been now reprinted in this same FCAA issue, under the 
kind permission of Birkhauser Verlag AG. 
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Figure 1: Plots of the function with /x = 1/4, 1/2,3/4, 1 versus t; top: 
hnear scales (0 < t < 5); bottom: logarithmic scales (10~^ <t< 10^). 

The plots of the functions '^{t) and $(/!:) are shown in Figs. 1 and 2, respec- 
tively, for some rational values of the parameter yU, by adopting linear and 
logarithmic scales. 

For the use of the Hilfer intermediate derivative in fractional relaxation we 
refer to Hilfer himself, see [67], p. 115, according to whom the Mittag-LefHer 
type function 

„(t)=t(i--)(M-i)i^^_^^^(^_^)(_tM), t>o, (2.10) 
is the solution of 

Df'^M(t) = t>0, with lim -^^^'-"^ = 1 . (2.11) 
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In fact, Hilfer has shown that the Laplace transform of the solution of (2.11) 
is 

s^^ + 1 

so, as a consequence of (2.7), we find (2.10). For plots of the function in 
(2.10) we refer to [115j. 
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3 Constitutive equations of fractional order 
in viscoelasticity 

In this section we present the fundamentals of hnear Viscoelasticity restrict- 
ing our attention to the one-axial case and assuming that the viscoelastic 
body is quiescent for all times prior to some starting instant that we assume 
as t — 0. 

For the sake of convenience both stress a{t) and strain e(i) are intended to 
be normalized, i.e. scaled with respect to a suitable reference state {a"o , eo} . 
After this necessary introduction we shall consider the main topic concerning 
viscoelastic models based on differential constitutive equations of fractional 
order. 

3.1 Generalities 

According to the linear theory, the viscoelastic body can be considered as a 
linear system with the stress (or strain) as the excitation function (input) and 
the strain (or stress) as the response function (output). In this respect, the 
response functions to an excitation expressed by the Heaviside step function 
0(t) are known to play a fundamental role both from a mathematical and 
physical point of view. We denote by J{t) the strain response to the unit 
step of stress, according to the creep test and by G{t) the stress response to 
a unit step of strain, according to the relaxation test. 

The functions J{t) , G{t) are usually referred to as the creep compliance and 
relaxation modulus respectively, or, simply, the m,aterial functions of the vis- 
coelastic body. In view of the causality requirement, both functions are 
causal, i.e. vanishing for t < 0. 

The limiting values of the material functions for t — > 0+ and t — > -|-oo 
are related to the instantaneous (or glass) and equilibrium behaviours of 

the viscoelastic body, respectively. As a consequence, it is usual to denote 
Jg := J{0^) the glass compliance, J,. := J(+oo) the equilibrium compliance, 
and Gg := G{0^) the glass modulus Ge '.= ^(+00) the equilibrium modulus. 
As a matter of fact, both the material functions are non-negative. Further- 
more, for < i < -|-oo , J{t) is a nan decreasing function and G{t) is a non 
increasing function. 

The monotonicity properties of J{t) and G{t) are related respectively to 
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the physical phenomena of strain creep and stress relaxatioi^. Under the 
hypotheses of causal histories, we get the stress-strain relationships 



e(t) 



r J(t - r) da{T) = a(0+) J(t) + / J(t - r) ^a{T) dr , 
Jo- Jo dr 

a(t) = r G(t - r) deir) = e(0+) G(t) + fcit - r) ^e(r) dr , 
JO- JO dr 



(3.1) 



where the passage to the RHS is justified if differentiability is assumed for 
the stress-strain histories, see also the excellent book by Pipkin j^Sj. Be- 
ing of convolution type, equations (3.1) can be conveniently treated by the 
technique of Laplace transforms so they read in the Laplace domain 



e(s) = s J{s) a{s) , cr(s) = s G{s) e(s) 
from which we derive the reciprocity relation 

s J{s) ^ 



sG{s) 



(3.2) 



(3.3) 



Because of the limiting theorems for the Laplace transform, we deduce that 
Jg = 1/Gg, Je = l/Ge, with the convention that and +oo are reciprocal 
to each other. The above remarkable relations allow us to classify the vis- 
coelastic bodies according to their instantaneous and equilibrium responses 
in four types as stated by Caputo & Mainardi in their 1971 review paper [29] 
in Table 3.1. 



Type 


J9 


Je 


Gg 


Ge 


1 


> 


< 00 


< 00 


> 


11 


> 


= 00 


< 00 


= 


111 


= 


< 00 


= 00 


> 


IV 


= 


= 00 


= 00 


= 



Table 3.1 The four types of viscoelasticity. 



■^For the mathematical conditions that the material functions must satisfy to agree 
with the most common experimental observations (physical realizability) we refer to the 
recent paper by Hanyga [63] and references therein. We also note that in some cases the 
material functions can contain terms represented by generalized functions (distributions) 
in the sense of Gel'fand-Shilov [49] or pseudo-functions in the sense of Doetsch [36] 
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We note that the viscoelastic bodies of type I exhibit both instantaneous 
and equiUbrium elasticity, so their behaviour appears close to the purely 
elastic one for sufficiently short and long times. The bodies of type II and 
IV exhibit a complete stress relaxation (at constant strain) since Ge = and 
an infinite strain creep (at constant stress) since Jg = oo , so they do not 
present equilibrium elasticity. Finally, the bodies of type III and IV do not 
present instantaneous elasticity since Jg — (Gg — oo). Other properties 
will be pointed out later on. 

3.2 The mechanical models 

To get some feeling for linear viscoelastic behaviour, it is useful to consider 
the simpler behaviour of analog mechanical models. They are constructed 
from linear springs and dashpots, disposed singly and in branches of two 
(in series or in parallel). As analog of stress and strain, we use the total 
extending force and the total extension. We note that when two elements 
are combined in series [in parallel], their compliances [moduli] are additive. 
This can be stated as a combination rule: creep compliances add in series, 
while relaxation moduli add in parallel. 

The mechanical models play an important role in the literature which is justi- 
fied by the historical development. In fact, the early theories were established 
with the aid of these models, which are still helpful to visualize properties and 
laws of the general theory, using the combination rule. Now, it is worthwhile 
to consider the simple models of Fig. 3 providing their governing stress-strain 
relations along with the related material functions. 

The spring (Fig. 3a) is the elastic (or storage) element, as for it the force is 
proportional to the extension; it represents a perfect elastic body obeying the 
Hooke law (ideal solid). This model is thus referred to as the Hooke model. 
If we denote by m the pertinent elastic modulus we have 



In this case we have no creep and no relaxation so the creep compliance 
and the relaxation modulus are constant functions: J{t) = Jg = Jg — 1/m; 



The dashpot (Fig. 3b) is the viscous (or dissipative) element, the force 
being proportional to rate of extension; it represents a perfectly viscous body 



Hooke model : a{t)=me{t), and 




(3.4) 



G{t) = Gg = Ge^ 1/m. 
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obeying the Newton law (perfect liquid). This model is thus referred to as 
the Newton model. If we denote by h the pertinent viscosity coefficient, we 
have 



Newton model : a{t) = bi — , and < ^ ^ 5^ ' (3.5) 

[G{t) = h6it). 

In this case we have a linear creep J{t) — J^t and instantaneous relaxation 
G{t) = 6{t) with G. = 1/ J+ = bi. 

We note that the Hooke and Newton models represent the limiting cases of 
viscoelastic bodies of type / and IV, respectively. 




Figure 3: The representations of the basic mechanical models: a) spring for 
Hooke, b) dashpot for Newton, c) spring and dashpot in parallel for Voigt, 
d) spring and dashpot in series for Maxwell. 



A branch constituted by a spring in parallel with a dashpot is known as the 
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Voigt model (Fig. 3c). We have 



Voigt model : a{t) = me{t) + 6] 



de 



(3.6) 



and 



J{t) = Ji 



-tjr. 



Jl = — , = — 

m m 



G{t) = Ge + G-5{t), Ge = m,G- 
where is referred to as the retardation time. 

A branch constituted by a spring in series with a dashpot is known as the 
Maxwell model (Fig. 3d). We have 



Maxwell model 



da de 
a(t) + ai- = 6,- 



(3.7) 



and 



J{t) = Jg + J+t, Jg 



G{t) =Gie 



Gi 



ai 



J4 



1 

b[ 

ai 



where To- is is referred to as the the relaxation time. 

The Voigt and the Maxwell models are thus the simplest viscoelastic bodies 
of type /// and //, respectively. The Voigt model exhibits an exponential 
(reversible) strain creep but no stress relaxation; it is also referred to as the 
retardation element. The Maxwe// model exhibits an exponential (reversible) 
stress relaxation and a linear (non reversible) strain creep; it is also referred 
to as the relaxation element. 

Based on the combination rule, we can continue the previous procedure in 
order to construct the simplest models of type / and IV that require three 
parameters. 

The simplest viscoelastic body of type I is obtained by adding a spring either 
in series to a Voigt model or in parallel to a Maxwell model (Fig. 4a and 
Fig 4b, respectively). So doing, according to the combination rule, we add a 
positive constant both to the Voigt-like creep compliance and to the Maxwell- 
like relaxation modulus so that we obtain Jg > and Ge > . Such a model 
was introduced by Zener [120J with the denomination of Standard Linear 
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Solid {S.L.S.). We have 



Zener model 



d' 



a{t) 



e{t), 



(3.8) 



and 



J{t) ^Jg + Jl 



G{t) ^Ge + Gi e" 



-t/Te 



T T ---Z± -h 

9 — L ' 1 ; ' ' 



-t/r. 



m 
bi 



bi 



m 



m, G\ — m, Ta — a\ 

ai 



We point out the condition < m < 6i/ai in order Ji,Gi be positive and 
hence < Jg < Jg < oo and < G^ < Gg < oo . As a consequence, we note 
that, for the S.L.S. model, the retardation time must be greater than the 
relaxation time, i.e. < To- < < oo . 




a) b) c) d) 

Figure 4: The mechanical representations of the Zener [a), b)] and anti- 
Zener[c), d)] models: a) spring in series with Voigt, b) spring in parallel 
with Maxwell, c) dashpot in series with Voigt, d) dashpot in parallel with 
Maxwell. 

Also the simplest viscoclastic body of type IV requires three parameters, 
i. e. oi , &i , 62 ; it is obtained adding a dashpot either in series to a Voigt 
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model or in parallel to a Maxwell model (Fig. 4c and Fig 4d, respectively). 
According to the combination rule, we add a linear term to the Voigt-like 
creep compliance and a delta impulsive term to the Maxwell-like relaxation 
modulus so that we obtain Jg = oo and Gg = oo . We may refer to this model 
to as the anti-Zener model. We have 



anti — Zener model : 



d' 



a{t) 



d d^' 
^'dt+^'d^ 



e{t) , (3.9) 



and 



J{t) = J+t + Ji 



G{t) = G_S{t) + Gie 



1-1. T - h —bi 

bi bi bf bi 

— — , Cji — n, To- — tti- 

ai ai ai 



We point out the condition < 62/^'i < Qi in order Ji,Gi be positive. As 

a consequence, we note that, for the anti-Zener model, the relaxation time 
must be greater than the retardation time, i.e. < < To- < oo , on the 
contrary of the Zener (S.L.S.) model. 

In Fig. 4 we exhibit the mechanical representations of the Zener model (3.8), 
see a), b), and of the anti-Zener model (3.9)), see c), d). 

Based on the combination rule, we can construct models whose material 
functions are of the following type 



J{t) = Jg + Y.Jn 



+ J+t, 



G{t) ^Ge + J2Gne -^/^-.- + S{t) , 



(3.10) 



where all the coefficient are non-negative, and interrelated because of the 
reciprocity relation (3.3) in the Laplace domain. We note that the four types 
of viscoelasticity of Table 3.1 are obtained from Eqs. (3.10) by taking into 
account that 



' Jp < oo 



Gg < oo 



J+ = , Je = oo <^=r- J+ 7^ , 
G'_ = 0, Gg^oo <^ G^^O. 



(3.11) 
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Appealing to the theory of Laplace transforms, we write 



sJ{s) = J, + Y.- 

'71 -L 



Jn 



+ S Te .„ S 
n T '(7,1 



(3.12) 



G^s 



where we have put /? = J2n Gn ■ 

Furthermore, as a consequence of (3.12), sJ{s) and sG{s) turn out to be 
rational functions in C with simple poles and zeros on the negative real axis 
and, possibly, with a simple pole or with a simple zero at s = , respectively. 



In these cases the integral constitutive equations (3.1) can be written in 
differential form. Following Bland |i)J with our notations, we obtain for these 
models 



ait) 



m 



k=l 



(3.13) 



where q and p are integers with q = p or q = p+ l and m, ak, hk are non- 
negative constants, subjected to proper restrictions in order to meet the 
physical requirements of realizability. The general Eq. (3.13) is referred to 
as the operator equation of the mechanical models. 

In the Laplace domain, we thus get 



sJ{s) 



Pis) 



sG{s) Q{s) 



where 



P{s) 
Q{s) 



1 + £ afc s\ 
fc=i 



k=l 



(3.14) 



with m > and q = p or q = p+ l. The polynomials at the numerator 
and denominator turn out to be Hurwitz polynomials (since they have no 
zeros for Re {s} > 0) whose zeros are alternating on the negative real axis 
(s < 0). The least zero in absolute magnitude is a zero of Q{s). The four 
types of viscoelasticity then correspond to whether the least zero is ( J+ 7^ 0) 
or is not ( J+ = 0) equal to zero and to whether the greatest zero in absolute 
magnitude is a zero of P{s) {Jg 7^ 0) or a zero of Q{s) {Jg = 0). 
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In Table 3.2 we summarize the four cases, which are expected to occur in the 
operator equation (3.13), corresponding to the four types of viscoelasticity. 



Type 


Order 


m 


■J9 




■J+ 




I 




> 




m 








II 


q = p 


= 


dp /hp 





l/bi 





III 


q = p + 1 


> 





m 





hq/ap 


IV 


q = p + 1 


= 








l/hi 


hq/ap 



Table 3.2: The four cases of the operator equation. 



We recognize that for p = 1, Eq. (3.13) includes the operator equations for 
the classical models with two parameters: Voigt and Maxwell, illustrated in 
Fig. 3, and with three parameters: Zener and anti-Zener, illustrated in Fig. 
4. In fact we recover the Voigt model (type III )for m > and p = 0,q = 1, 
the Maxwell model (type II) for m = and p — q — 1, the Zener model (type 
I) for m > and p — q — 1, and the anti-Zener model (type IV) for m = 
and p = l,q = 2. 

Remark. We note that the initial conditions at t = 0+, a^'^^O'^) with h = 
0, 1, . . .p — 1 and £^'^•'(0"'") with k = 0, 1, ... g — 1, do not appear in the 
operator equation but they are required to be compatible with the integral 
equations (3.1). In fact, since Eqs (3.1) do not contain the initial conditions, 
some compatibility conditions at t = O"*" must be implicitly required both for 
stress and strain. In other words, the equivalence between the integral Eqs. 
(3.1) and the differential operator Eq. (3.13) imphes that when we apply 
the Laplace transform to both sides of Eq. (3.13) the contributions from 
the initial conditions are vanishing or cancel in pair-balance. This can be 
easily checked for the simplest classical models described by Eqs (3.6)-(3.9). 
It turns out that the Laplace transform of the corresponding constitutive 
equations does not contain any initial conditions: they are all hidden being 
zero or balanced between the RHS and LHS of the transformed equation. 
As simple examples let us consider the Voigt model for which p — 0, q — 1 
and m > 0, see Eq. (3.6), and the Maxwell model for which p — q — 1 and 
m = 0, see Eq. (3.7). 

For the Voigt model we get sa{s) = mt{s) + h]se{s) — e(0"'")], so, for any 
causal stress and strain histories, it would be 

sJ(s) = — L— e(0+) = 0. (3.15) 
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We note that the condition e(0''") = is surely satisfied for any reasonable 
stress history since Jg = 0, but is not valid for any reasonable strain history; 
in fact, if we consider the relaxation test for which e{t) = 0(t) we have 
e(0+) = 1. This fact may be understood recalling that for the Voigt model 
we have Jg — and Gg — oo (due to the delta contribution in the relaxation 
modulus) . 

For the Maxwell model we get a{s) + a [sa{s) - cr(0+)] = h [se{s) - e(0+)], 
so, for any causal stress and strain histories it would be 

sJ(s) = T + T- acx(0+) = 6e(0+) . (3.16) 

OS 

We now note that the condition aa{0'^) — 66(0+) is surely satisfied for any 
causal history both in stress and in strain. This fact may be understood 
recalling that for the Maxwell model we have Jg > and Gg = 1/ Jg > 0. 

Then we can generalize the above considerations stating that the compat- 
ibility relations of the initial conditions are valid for all the four types of 
viscoelasticity, as far as the creep representation is considered. When the 
relaxation representation is considered, caution is required for the types III 
and IV, for which, for correctness, we would use the generalized theory of 
integral transforms suitable just for dealing with generalized functions. 

3.3 The time spectral functions 

From the previous analysis of the classical mechanical models in terms of a 
finite number of basic elements, one is led to consider two discrete distribu- 
tions of characteristic times (the retardation and the relaxation times), as 
stated in (3.10). However, in more general cases, it is natural to presume 
the presence of continuous distributions, so that, for a viscoelastic body, the 
material functions turn out to be of the following form: 

J{t) = Jg + a /o^ i?e(r) (l - e-t/A dr + J+t, 

^ ^ (3.17) 

^ G{t) = Ge + /3 /o°° Ra{r) e-^/^ dr + G_ 5{t) . 

where all the coefficients and functions are non- negative. The function Re{T) 
is referred to as the retardation spectrum while Ra{T) as the relaxation spec- 
trum. For the sake of convenience we shall replace the suffix e or r with * to 
denote anyone of the two spectra that we refer simply to as the time-spectral 
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function. We require R^, (r) to be locally integrable in R"*" , with the supple- 
mentary normalization condition /q°° R{t) dr = 1 ii the integral in R"*" turns 
out to be convergent. 

The discrete distributions of the classical mechanical models, see (3.10), can 
be easily recovered from (3.17); in fact, assuming a 7^ , /5 7^ , we have to 
put 

" " " (3.18) 

We devote particular attention to the time-dependent contributions to the 
material functions (3.17) which are provided by the continuous spectra, i.e. 

qf(t) ■.= a r Reir) (l-e^^/A dr , 

^ ^ (3.19) 

$(t) ■= p / i?^(r)e"*/^c/r. 
Jo 

We recognize that \l/(t) (that we refer as the creep function with spectrum) is a 
non- decreasing, non-negative function in R"*" with limiting values \l/(0"^) = 0, 
\E'(+oo) = a or 00, whereas $(t) (that we refer as the relaxation function 
with spectrum) is a non-increasing, non-negative function in R"*" with limiting 
values $(0"*") = (3 or 00, $(-|-oo) = 0. More precisely, in view of their spectral 
representations (3.19), we have 

^(t)>o, (-ir^<o, ^ ^ 

t>0, n = l,2,.... (3.20) 

In other words, $(t) is a completely monotonic function and \l/(t) is a Bern- 
stein function (namely a non-negative function with a completely monotonic 
derivative). These properties have been investigated by several authors, in- 
cluding Molinari [87] and more recently by Hanyga [63]. The determina- 
tion of the time-spectral functions starting from the knowledge of the creep 
and relaxation functions is a problem which can be formally solved through 
the Titchmarsh inversion formula of the Laplace transform theory, see e.g. 
[621, [29]. 
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3.4 Fractional viscoelastic models 



The straightforward way to introduce fractional derivatives in hnear vis- 
coelasticity is to replace the first derivative in the constitutive equation (3.5) 
of the Newton model with a fractional derivative of order v G (0, 1), that, 
being e(0"'") = 0, may be intended both in the Riemann-Liouville or Caputo 
sense. Some people call the fractional model of the Newtonian dashpot with 
the suggestive name 'pot. we prefer to refer such model to as Scott-Blair 
model, to give honour to the scientist who already in the middle of the past 
century proposed such a constitutive equation to explain a material prop- 
erty that is intermediate between the elastic modulus (Hooke solid) and the 
coefficient of viscosity (Newton fluid), see e.g. [1111 11131 [114j . We note that 
Scott-Blair was surely a pioneer of the fractional calculus even if he did not 
provide a mathematical theory accepted by mathematicians of his time! 

The use of fractional calculus in linear viscoelasticity leads to generalizations 
of the classical mechanical models: the basic Newton element is substituted 
by the more general Scott-Blair element (of order z/). In fact, we can construct 
the class of these generalized models from Hooke and Scott-Blair elements, 
disposed singly and in branches of two (in series or in parallel) . Then, extend- 
ing the procedures of the classical mechanical models (based on springs and 
dashpots), we will get the fractional operator equation (that is an operator 
equation with fractional derivatives) in the form which properly generalizes 
(3.13), I.e. 



k=l 



ait) 



1 



k=l 



df" 



e{t) , Uk = k + u-l. (3.21) 



so, as a generalization of (3.10), 



At) = A + i: Jn {1 - K [-{tKuY]} + J+YTI 



G{t) = a + E G'n [-{tKnY] + 



(3.22) 



rri - u) ' 



where all the coefficient are non-negative. Of course, for the fractional oper- 
ator equation (3.21) the four cases summarized in Table 3.2 are expected to 
occur in analogy with the operator equation (3.13). 
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Figure 5: Spectral function R^{t) for z/ = 0.25, 0.50, 0.75, 0.90 and = 1. 



We conclude with some considerations on the presence of the Mittag-Leffler 
function in the material functions of the fractional models. For this pur- 
pose let us consider the following creep and relaxation functions with the 
corresponding time-spectral functions 



^(t) = a {1 -E,[-{t/T,y]} = a J\{t) (^1-e-^/^) dr , 

m = PE,[-it/T^y]=P rR„{T)e-^l^ dr. 

Jo 



(3.23) 



Following Caputo and Mainardi [281 EH] and denoting with star the suffixes 
e , 0" , we obtain 



RJr) 



sm Z/TT 



71 T {t / + {r / T^) ^ + 2 COS z/vr 



(3.24) 



From Fig 5 reporting the plots of -R*(t) for some values of z/ we can easily 
recognize the effect of the variation of u on the character of the spectrum. For 
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< u < uq, where z/q ~ 0.736 is the solution of the equation u = sin uir, the 
spectrum -R*(r) is a decreasing function of r; subsequently, with increasing 
z/, it first exhibits a minimum and then a maximum; for z/ — »• 1 it becomes 
steeper and steeper near its maximum approaching a delta function. In fact 
we have lim i?*(r) = 6(t — r*), where we have set r* = 1, being r* the single 

retardation/relaxation time exhibited by the corresponding creep/relaxation 
function. Formula (3.24) was formerly obtained by Gross [BT] in 1947, when, 
in the attempt to eliminate the faults which a power law shows for the creep 
function, he proposed the Mittag-LefHer function as an empirical law for 
both the creep and relaxation functions. In their 1971 papers, Caputo and 
Mainardi derived the same result by introducing into the stress-strain rela- 
tions a Caputo derivative that implies a memory mechanism by means of a 
convolution between the first-order derivative and a power of time, see (1.6'). 



In fractional viscoelasticity governed by the operator equation (3.21) the 
corresponding material functions are obtained by using the combination rule 
valid for the classical mechanical models. Their determination is made easy if 
we take into account the following correspondence principle between the clas- 
sical and fractional mechanical models, as introduced by Caputo & Mainardi 
in 1971. Taking < u < 1 such correspondence principle can be formally 
stated by the following three equations where transitions from Laplace trans- 
form pairs are outlined: 

1 r 1 , , 

where r > and E,y denotes the Mittag-LefHer function of order u. 

Remark. We note that the initial conditions at t = O"*" for the stress and 
strain do not explicitly enter into the fractional operator equation (3.21) if 
they are taken in the same way as for the classical mechanical models re- 
viewed in Subsection 3.2. This means that the approach with the Caputo 
derivative, which requires in the Laplace domain the same initial conditions 
as the classical models is quite correct. However, if we assume the same 
initial conditions, the approach with the Riemann-Liouville derivative pro- 
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vides the same results since, in view of the corresponding Laplace transform 
rule (1.15), the initial conditions do not appear in the Laplace domain. The 
equivalence of the two approaches has also been noted for the fractional Zener 
model in a recent note by Bagley [5|. We refer the reader to the paper by 
Heymans and Podubny [65] for the physical interpretation of initial condi- 
tions for fractional differential equations with Riemann-Liouville derivatives, 
especially in viscoelasticity. In such field, however, we prefer to adopt the 
Caputo derivative since it requires the same initial conditions as in the clas- 
sical pointed out in [81]: by the way, for a physical view point, these 
initial conditions are more accessible than those required in the more general 
Riemann-Liouville approach, see (1-14). 

4 Some historical notes 

In this section we provide some historical notes on how the Caputo derivative 
came out as a contrast to the Riemann-Liouville derivative and a sketch about 
the generic use of fractional calculus in viscoelasticity. 

4.1 The origins of the Caputo derivative 

Here we find it worthwhile and interesting to say something about the com- 
monly used attributes to Riemann and Liouville and to Caputo for the two 
types of fractional derivatives, that have been discussed. 

Usually names are given to pay honour to some scientists who have provided 
main contributions, but not necessarily to those who have as the first in- 
troduced the corresponding notions. Surely Liouville, e.g. [75l [76] (starting 
from 1832), and then Riemann [105] (as a student in 1847) have given im- 
portant contributions towards fractional integration and differentiation, but 
these notions had previously a story. As a matter of fact it was Abel who 
solved his celebrated integral equations by a fractional integration of order 
1/2 and a G (0, 1), respectively in 1823 and 1826, see [HIS]. So, Abel, using 
the operators that nowadays are ascribed to Riemann and Liouville, preceded 
these eminent mathematicians by at least 10 years! 

We remind that the Laplace transform rule (1.13) was practically the start- 
ing point of Caputo [T2l [T3] in defining his generalized derivative in the 
late Sixties of the past century, and ignoring the existence of the classical 
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Riemann-Liouville derivative. Please keep in mind that the first treatise de- 
voted to the so-called fractional calculus appeared only in 1974, published 
by Oldham & Spanier [91] who were unaware of the alternative form (1.6) 
of the fractional derivative and of its property (1.13) with respect to the 
Laplace transform. However the form used by Caputo is found in a paper by 
Liouville himself as recently noted by Butzer and Westphal [TU] but Liouville 
disregarded this notion because he did not recognize its role. 

As far as we know, up to the middle of the past century the authors did 
not take care of the difference between the two forms (1.5)- (1.6) and of the 
possible use of the alternative form (1.6). Indeed, in the classical book on 
Differential and Integral Calculus by the eminent mathematician R. Courant 
the two forms of the fractional derivative were considered as equivalent, see 
[3l] . pp. 339-341. Only in the late Sixties it seems that the relevance of 
the alternative form was recognized. In fact, in 1968 Dzherbashyan and 
Nersesian [38j used the alternative form for dealing with Cauchy problems 
of differential equations of fractional order. In 1967, just one year earlier, 
Caputo [12], see also [13], used this form to generalize the usual rule for 
the Laplace transform of a derivative of integer order and to solve some 
problems in Seismology. Soon later, this derivative was adopted by Caputo 
and Mainardi in the framework of the theory of Linear Viscoelasticity, see 

[2HIES]. 

Starting with the Seventies many authors, very often ignoring the works 
by Dzherbashian-Nersesian and Caputo-Mainardi, have re-discovered and 
used the alternative form, recognizing its major utility for solving physi- 
cal problems with standard (namely integer order) initial conditions. Al- 
though there appeared several papers by different authors, including Ca- 
puto [ll[T5l[TSl[I71[l8l[T9l[20l[2ll[22l[2a Gorenflo et al. 
[551ESIEZIEH], Kochubei [HI El, Mainardi [ISl ESI EOl ES] , where the alter- 
native derivative was adopted, it was mainly with the 1999 book by Podlubny 
[96] that it became popular: indeed, in that book it was named as Caputo 
derivative. 

The notation adopted here was introduced in a systematic way by us in our 
1996 CISM lectures [57], partly based on the 1991 book on Abel Integral 
Equations by Gorenflo & Vessella [59] . 
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4.2 Fractional calculus in viscoelasticity in the XX-th 
century 

Starting from the past century a number of authors have (imphcitly or 
exphcitly) used the fractional calculus as an empirical method of describ- 
ing the properties of viscoelastic materials: Gemant (SOj [52], Scott-Blair 
|114[ I113t I112j . Gerasimov [53] were early contributors. Particular mention 
is due to the theory of hereditary solid mechanics developed by Rabotnov 
|1U2[ \WS\ llU4j . see also |1U7] . that {implicitly) requires fractional derivatives. 



In the late Sixties, Caputo [TTl [121 [IS] and then Caputo and Mainardi [28l 
suggested that derivatives of fractional order (of Caputo type) could be 
successfully used to model the dissipation in seismology and in metallurgy. 

From then, applications of fractional calculus in viscoelsticity were consid- 
ered by an increasing number of authors. Restricting our attention to a few 
significant papers published in the past century, let us quote the contribu- 
tions by Bagley & Torvik H El HII] , Caputo [H [161 [HI [201 [251 [27] , Friedrich 
& associates HQl [HI [H [l3l [B], Graffi [60], Gaul, Kempfle & associates 
[HI [13 [IHl [17], Heymans [M], Koeller [73l [H], Mainardi [THl [80], Meshkov 
and Rossikhin [M] , Nonnenmacher & associates [SH [121 [Si] , Pritz [SSI [TUP] , 
Rossikhin and Shitikova [106j . 

Additional references up to nowadays can be found in the huge (even not 
exhaustive) bibliography of the forthcoming book by Mainardi [H2] • 

We like to recall the formal analogy between the relaxation phenomena in 
viscoelastic and dielectric bodies; in this respect the pioneering works by Cole 
and Cole [221 133] in 1940's on dielectrics can be considered as precursors of 
the implicit use of fractional calculus in that area, see e.g. [86] . 



Appendix: The functions of Mittag-Leffler type 

A.l. The classical Mittag-Leffler function 

The Mittag-Leffler function E^{z) (with /i > 0) is an entire transcendental 
function of order l//i, defined in the complex plane by the power series 

oo k 
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It was introduced and studied by the Swedish mathematician Mittag-Leffler 
at the beginning of the XX-th century to provide a noteworthy example 
of entire function that generahzes the exponential (to which it reduces for 

Details on this function can be found e.g. in the treatises by Davis [35] . 
Dzherbashyan [37|, Erdelyi et al. [SH], Kilbas et al. [02], Kiryakova [70] . 
Podlubny [HS], Samko et al. |lU8j . Sansone and Gerretsen [109j . Concerning 
earlier applications of the Mittag-Leffler function in physics let us quote the 
contributions by K. S. Cole, see [31], mentioned in the 1936 book by Davis, 
see [35], p. 287, in connection with nerve conduction, and by Gross, see [6T] . 
in connection with creep and relaxation in viscoelastic media. 

We note that the function E^{—x) {x > 0) is a completely monotonic function 
of X if < ;U < 1, as was formerly conjectured by Feller on probabilistic 
arguments and later in 1948 proved by Pollard [98] . This property still holds 
with respect to the variable t if we replace x hj Xf^ {t > 0) where A is 
a positive constant. Thus in its dependence on t the function At^) 
preserves the complete monotonicity of the exponential exp(— At): indeed, 
for < yU < 1 it is represented in terms of a real Laplace transform (of a 
real parameter r) of a non-negative function (that we refer to as the spectral 
function) 

^ / , „x 1 /■°° —rf Ar'^"^ sin(u7r) , , ^ , 

EJ-Xtn = - e ^ )^ -dr. (A.2) 

'^^ ' Jo A2 + 2Ar/^ cos(;U7r) +r2M ^ ' 

We note that as /i — 1~ the spectral function tends to the generalized Dirac 
function 5{r — A). 

We point out that the Mittag-Leffler function (A.2) starts at t = as a 
stretched exponential and decreases for t — > cxd like a power with exponent 



tf" f At^ 1 
1-A— r~exp<^-— -) , t^O" 



r(i + /x) r(i + /.) 



(A3) 



t — » oo . 



I Ar(i-/i)' 

The integral representation (A.2) and the asymptotic (A. 3) can also be de- 
rived from the Laplace transform pair 



C{E,{-Xn,s] = -^. (A4) 
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In fact it it sufficient to apply the Titchmarsli theorem for contour integration 
in the complex plane (s = re^^) for deriving (A. 2) and the Tauberian theory 
{s ^ oo and s — 0) for deriving (A. 3). 
If /X = 1/2 we have for t > 0: 

Ei/2{-\Vt) = e erfc(Av^) ~ l/(Av^) , t ^ oo , {A.5) 
where erfc denotes the complementary error function, see e.g. [3]. 



A. 2. The generalized Mittag-Leffler function 

The Mittag-Leffler function in two parameters E^^v{z) (3?{/i} > 0, z/ G C) is 
defined by the power series 

oo k 

It generalizes the classical Mittag-Leffler function to which it reduces for 
V = 1. It is an entire transcendental function of order on which 

the reader can inform himself by again consulting the treatises cited for the 
classical Mittag-Leffler function. 

The function E^y{—x) {x > 0) is completely monotonic in x if < /x < 1 
and > /i, see e.g. [89 l [90 t [TTO] . Again this property still holds with respect 
to t if we replace the variable x by At'^ where A is a positive constant. In 
this case the asymptotic representations a.s t ^ and t — » +oo read 



A— r, ^^0^ 



T{u) r(z/ + /i) 
1 t-^'+^-^ 



{A.7) 



oo. 



A r(z/-/i) 

We point out the Laplace transform pair, see [96] , 

C{e-^E,A-Xn,s} = ^^, (A8) 

+ A 

with /i, z/ G R"*". By aid of this Laplace transform, with 0</i = z/<l, we 
can obtain the useful identity 

r(i-'^)E^,^(-AtO = -^^E^(-At^) , </.<!. (A9) 
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To see this it is sufficient to write 



1 



1 r si^-^ 



- 1 



(^.10) 



X s^' + X 



and invert tlie Laplace transforms. Of course tlie identity (A. 9) can be proved 
directly by differentiating term by term the power series of the classical 
Mittag-Leffler function, but, as often in matters of fractional calculus, it 
is simpler to work with the Laplace transform technique. 
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